use dataset, clear

net install grc1leg, replace from("http://www.stata.com/users/vwiggins/")

* Collapse

collapse ft_dist polar_dif, by(year)
sum year
gen time = (year-r(min))/(r(max)-r(min))
orthpoly time, generate(time1 time2) degree(2)

* Mean difference

twoway ///
(scatter ft_dist year) ///
(fpfit ft_dist year) ///
, ///
title("Mean difference") ///
xtitle("Year") ///
ytitle("Average absolute difference") ///
ylab(, angle(horiz) nogrid) ///
legend(rows(1) order(2 1) label(2 "Estimates") label(1 "Polynomial trend")) ///
scheme(s2mono) ///
graphregion(fcolor(white) lcolor(white)) ///
name(g1, replace) ///
nodraw

* Bimodality

twoway ///
(scatter polar_dif year) ///
(fpfit polar_dif year) ///
, ///
title("Bimodality") ///
xtitle("Year") ///
ytitle("Esteban--Ray polarization index") ///
ylab(, angle(horiz) nogrid) ///
legend(rows(1) order(2 1) label(2 "Estimates") label(1 "Polynomial trend")) ///
scheme(s2mono) ///
graphregion(fcolor(white) lcolor(white)) ///
name(g2, replace) ///
nodraw

* Combine

grc1leg g1 g2, ///
rows(1) ///
scheme(s2mono) ///
graphregion(fcolor(white)) ///
name(g12, replace)

* Export
gr export "figure4.eps", replace
